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OV ABSTRACT 

O' 
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TeV 7 rays from distant astrophysical sources are attenuated due to electron- 
positron pair creation by interacting with ultraviolet /optical to infrared photons 

" . which fill the universe and are collectively known as the extra-galactic background 

light (EBL). We model the ~0.1-10 eV starlight component of the EBL derived 
from expressions for the stellar initial mass function, star formation history of 
the universe, and wavelength-dependent absorption of a large sample of galaxies 
q . in the local universe. These models are simultaneously fitted to the EBL data as 

well as to the data on the stellar luminosity density in our local universe. We find 
that the models with modified Salpeter A initial mass function together with Cole 
et al. (2001) or Hopkins & Beacom (2006) star formation history best represent 
available data. Since no dust emission is included, our calculated EBL models 
can be interpreted as the lower limits in the ~0.1-1 eV range. We present simple 

pq ! analytic fits to the best-fit EBL model evolving with redshift. We then proceed 

to calculate 7-ray opacities, and absorption of ~ 10-300 GeV 7-rays coming from 
different redshifts. We discuss implications of our results for the Fermi Gamma 
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Ray Space Telescope and ground-based Air Cherenkov Telescopes. 

O 



% 



Subject headings: stars: formation — stars: fundamental parameters — stars: lu- 
minosity function, mass function — dust, extinction — diffuse radiation — gamma 
rays: observations 



1. Introduction 

Stars are the dominant sources of electromagnetic radiation in the universe after the 
cosmic microwave background (see, e.g., Fukugita & Peebles 2004). They emit radiation 
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longward from ultraviolet to infrared wavelengths. However, photons with wavelength < 
2 /im are highly absorbed by the dust in the host galaxies and only a fraction of the radiation 
emitted by the stars escape to the inter-galactic medium and form a diffuse background or 
EBL (see, e.g., Baldry & Glazebrook 2003; Driver et al. 2008). The dust in the host galaxies, 
heated by the starlight, also radiate in the infrared wavelengths and contribute to the EBL 
density at ~ 10 12 Hz. It is the direct starlight component, < 2 /im or > 0.1 eV, that affects 
the propagation of < 5 TeV 7-rays from distant sources. Indeed, the very soft spectral 
energy distribution (dN/dE oc E~ r ) with T > 3 observed from several TeV blazars at high 
redshift (z > 0.1) such as PKS 2155-304 (Aharonian et al. 2005), H 2356-309 (Aharonian 
et al. 2006a); 1ES 1218+304 (Albert et al. 2006); 1ES 1101-232 (Aharonian et al. 2006b); 
0347-121 (Aharonian et al. 2007), 1ES 1011+496 (Albert et al. 2007) and 3C 279 (Albert et 
al. 2008), and their cutoff at > 1 TeV are hints that high energy 7 rays from these sources 
are absorbed by the EBL UV/optical photons (Persic & de Angelis 2008). Lower energy 
(<TeV) 7-rays from high redshift sources such as gamma-ray bursts (GRBs) and blazars 
can also probe the EBL starlight component. 

Calculation of the opacity of the universe to 7-rays by 77 — *■ e + e~ process dates back 
to Nishikov (1961), followed by Gould & Shreder (1966) and Fazio & Stecker (1970). More 
recently Malkan & Stecker (1998, 2001); Primack et al. (1999); Kneiske, Mannheim & 
Hartmann (2002); Kneiske et al. (2004); Primack, Bullock and Somerville (2005); Stecker, 
Malkan & Scully (2006) calculated EBL models adopting either a phenomenological approach 
or Monte Carlo galaxy formation code. These models trace the general trend of the data, 
which may be fitted with a combination of two or more modified blackbody spectra for its two 
distinct peaks at the infrared and optical wavebands (Dermer 2007). Significant uncertainty 
in data and large dispersion among models led to an indirect method to constrain the EBL, 
namely by estimating change in spectral slope from distant TeV blazars due to 77 absorption 
(Stecker & de Jager 1993; Stanev & Franceschini 1998; Mazin & Raue 2007). However, such 
a method generally does not include possible absorption at the source (see, e.g., Reimer 
2007) and presumes a source spectrum. 

In this paper, we build models of the EBL starlight component (~0.1-10 eV) directly 
from the stellar thermal surface radiation. Emission from an individual star during its main- 
sequence lifetime is well approximated as a blackbody with a mass-dependent temperature. 
The post-main-sequence lifetime of a star is very short compared to its main-sequence lifetime 
and their contribution at the UV-optical wave bands is not significant. They can, however, 
contribute significantly to longer wavelengths due to their increased luminosity in the post- 
main-sequence phase (Finke et al., in preparation). Only a small fraction of the stars with 
mass > 8M© produce supernovae and even a smaller fraction produce GRBs. Emission from 
these sources dominate the diffuse MeV background (Watanabe et al. 1999; Ruiz-Lapuente, 
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Casse & Vangioni-Flam 2001; see, however, Strigari et al. 2005; Inoue, Totani & Ueda 2008). 
Emission from quasars and AGNs, on the other hand, dominate the diffuse X-ray background 
(Mushotzky et al. 2000). Contributions from all these sources add only a small fraction to 
the total cosmic electromagnetic energy density in the ~ 0.1-10 eV range, and we also ignore 
that. The estimated lifetimes of individual stars depend on their masses and the assumed 
cosmology, which is the standard ACDM with (h, fi m , Qa) = (0.7, 0.3, 0.7), and the Hubble 
constant H = 70h /i km s _1 Mpc -1 . Summing over contributions from stars of all masses 
formed in the history of the universe then gives us the diffuse emission or EBL. 

A sum over contributions from individual stars radiating at a given redshift corresponds 
to the luminosity density or the stellar energy emissivity of the universe at that redshift. 
The initial mass function (IMF), which is the distribution of stars by mass, and the star 
formation rate (SFR), which is the mass that forms stars per unit comoving volume per unit 
time, are two uncertain but related parameters in our calculation. We form classes of models 
by choosing different combinations of these parameters and compare, in the UV-optical band, 
the luminosity density data of the local universe found from the surveys of nearby galaxies. 
The same models are then compared with EBL data. Note that there are no adjustable free 
parameters in our calculation once we choose a particular model. 

Finally we use one of our best-fit models to calculate the e ± pair production opacities in 
the ~10-300 GeV energy range at different redshifts. These results are applicable to high- 
energy emission from distant sources such as GRBs and blazars detected by the currently 
operating Fermi Gamma Ray Space Telescope and Air Cherenkov Telescopes such as HESS, 
MAGIC and VERITAS. 

In Sec. [2] we outline the formalism of our method and introduce different models in Sec. 
[3] upon which we base our EBL calculation. We report our results in Sec. HI compare with 
results from previous authors as well as calculate EBL evolution with redshift. We discuss 
implications of our results for 7-ray astronomy in Sec. [51 calculating 77 — > e + e~ absorption 
opacity (t 77 ) for high energy 7 rays from sources at different redshifts. Conclusions of the 
work for EBL and 7 ray absorption are given in Sec. [6j 



2. Formalism 

The differential number density (per unit volume) of thermal blackbody photons in the 
energy interval e to e + de at a given temperature T is 

dN 1 e 2 



dedV 7r 2 (hc) 3 exp(e/kT) 
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and the total number of photons emitted per unit energy and time intervals from a star of 
radius R is given by 

^ = rfe -£* . (2) 

The radius and temperature of a star depend on its mass M, and can be calculated in terms 
of the solar radius R Q , mass M and temperature T using various relations. We use a fit 
to the stellar mass-radius relation (Schmidt-Kaler 1982; Binney & Merrifield 1998) as 

f (M/M ) 08 ; O.1M < M < 1OM 
7 "" \ 1O 9 / 2O (M/M ) - 35 ; 1OM < M < 12OM . l ' 

The stellar mass-luminosity ratio is not well-known. To the first approximation, it is a 
single power law (SPL), given by 

L/L Q = (M/M ) 3 - 6 . (4) 

A more detailed relation is calculated by Bressan et al. (1993) and a simple broken power 
law (BPL) fit to that relation (Binney & Merrifield 1998) is given by 

r (M/M ) 4 - 8 ; M < 2M Q 
L/L Q = f L < 2 13 / 1O (M/M ) 3 - 5 ; 2M < M < 2OM (5) 

[ 2 2O1 / 5O 5 34 / 25 (M/M ) 2 - 14 ; M > 2OM . 

Although there are uncertainties, depending on the metallicity e.g., in the mass-to-light ratio, 
we assume /r, = 1 for our modeling purposes. Stars below 1M and above 2OM , in the 
Bressan et al. (1993) model, produce much less light than the model with a single power-law 
[equation (jl])] for the same masses. We derive a stellar mass-temperature relation from the 
luminosity L = ATrR 2 aT A , where a is the Stefan-Boltzmann constant. For the SPL relation 
in equation (J3J) between L — M, the temperature is 

r (M/M Q y/ 2 ; M < 1OM 
1 Q " \ lO- 9 / 4O (M/M ) - 725 ; M > 1OM , { ' 

following the break in equation Q. For the BPL relation of L — M in equation ([5]), the 
T — M relation is 

{(M/M ) - 8 ; O.1M < M < 2M 
2 13 / 4O (M/M ) - 475 ; 2M < M < 1OM 
2 i/io 5 -9/4o (m/Mq) o.7 . 10Mq < M < 20Mq U; 

2 39/5O 5 23/2OO (M//M0) O.36 . M > 2QM& _ 
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For reference, M Q = 1.99 x 10 33 g, L Q = 3.846 x 10 33 ergs s" 1 and T Q = 5777 K are solar 
mass, luminosity, and temperature respectively. 

One needs to take cosmology into account to calculate the total number of photons 
emitted from a star over cosmic time which was born at a past epoch. The relationship 
between the cosmic time and redshift is given by 

(dt/dz)- 1 = -h (i + z )^/n m (i + zf + n A . (8) 

The main sequence lifetime of a star with mass M is 

t**t Q (M/M & )/(L/L Q ), (9) 

where £ Q ~ 1 1 Gyr is the lifetime of the Sun. If a star of mass M was born at a redshift z 
then the redshift z&{M) at which it had evolved off the main sequence can be found from 
the inverse of the relation i*(M) = J (M) dz' \dt/dz'\ as 

z d (M,z) = -1 + (-(ft A /On) sech [(3/2)i7 £* \A^ 

+ tanh- 1 v/l + (n m /fi A )(l+^)3j J , (10) 

We have plotted this redshift in Fig. [1] for different M and z and setting it equal to zero 
if Zd{M) < 0. After calculating z& for a given (z,M) combination, we back-calculate z as 
a check. The results are identical. Note that high mass stars > 1OM evolve off the main 
sequence almost at the same redshift (z « 1 - 5) they were formed. On the contrary, a star 
with mass ~ 1M lives almost a Hubble time. Note that z d , following £*, depends on the 
L — M relation. For the curves in Fig. [1] we used both the SPL (dashed lines) and BPL 
(dotted lines) L — M relations given, respectively, by equations (J3J and (JSJ). 

Equation (Q represents an underlying source of uncertainty for high-mass stars and/or 
high redshift stars. Our naive estimate of stellar lifetime and luminosity are based on the 
model of the Sun, and assuming that stars are perfect blackbodies. Estimates of the Sun's 
age are based on the total amount of fusion material (4H + — * He ++ ) in the core. This 
amounts to about 10% of the Sun's mass. The estimates for high mass stars or stars with 
different metallicities are less precise. High mass stars often have winds which may reduce 
their net radiative output, however, they can also be approximated as perfect blackbodies 
because of their higly ionized surface. Since the fraction of high-mass stars is small, we 
neglect corrections to equation (jUJ) in this study. 

The integrated number of photons emitted by a star from its birth at redshift z to the 
present epoch, in the energy interval de = de'/(l + z'), can be calculated using equations (J2J), 
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©, © or ©, © or © and (JTOD as 



diV(e, M) 



dz' 



max{0, z d (M,z')} 



dt 
~dz' 



dN(e', M) 



de'dt 



1+z') 



(11) 



where the lower limit of the integration is set to zero if the star has survived to the present 
epoch. 

The number of stars formed at a redshift z depends on the initial mass function and 
star formation rate, both of which are important sources of uncertainty in our calculation 
and are discussed shortly. We assume a universal IMF which is normalized between mass 
(M min - M max ) = (0.1 - 12O)M as AT" 1 = JjJ*^* dM(dN/dM)M. The final number of 
photons reaching us from all stars created in the past, per unit energy interval, amounts to 
weighting equation ffTTT) by the normalized IMF, Af(dN/dM), and integrating over stellar 
mass (to take into account photons created by stars of all mass) and, after weighting by the 
SFR ij){z) in units of M yr _1 Mpc -3 , integrating over redshift z (to take into account stars 
created at all epochs). The spectral stellar radiation density is therefore given by 



dN(e, z = 0) 
de~dV 



M j dz" 

2=0 

dt 



dt 



x / dz' 

'max{0, z d (M,z')} 



dz' 



dz" 



M n 



ip{z") I dM 

,.diV(e',M) 



dN 



de'dt 



dM 
\+z'). 



(12) 



Here dN(e', M) /de'dt is given by equation (J2J) with the substituton e' = e(l + z'), and / eS c(e') 
is the escape fraction of photons from the host galaxy, considered further in the next section. 
Note that there is no free parameters in our model once we chose particular models of SFR, 
IMF and L - M relation. 



Equation f[T2"j) can be converted to the EBL energy density eu e (e.g., in units of ergs cm" 3 ) 
by multiplying with e 2 or to intensity el e (e.g., in units of W m~ 2 sr _1 ) by multiplying with 
e 2 c/47r. The differential photon number density in equation ffl2l or EBL energy density 
measured at present from a past epoch [z = Z\) can be transformed to the past epoch, i.e. 
comoving EBL density, with e\ = e(l + z\) and the comoving volume V\ = V/(l + zi) 3 as 



dNjeuzJ 
de\dV\ 

e x u ei 



(l + 2i 



2 dN(e, z = zi) 



dedV 



A2 dN(e,z = Zl ) 

^ 1 + Zl)e — dldv — 



(13) 



Note that we use z = z\ as the lower limit of the outer integration over redshift in equa- 
tion ffl2l to calculate dN(e, z = z^/dedV. 
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2.1. Energy output by the local universe 

The cosmic energy output or luminosity density eL e (e.g. in units of W Mpc -3 ) in 
starlight is generally found from the galaxy counts in the local universe (see, e.g., Baldry & 
Glazebrook 2003; Driver et al. 2008). This is equivalent to summing over energy output by 
individual stars which were born at a past epoch but still radiating today (z — Z\ « 0), and 
we calculate the comoving luminosity density by modifying equation (TT2|) as 



/■M max / AfJ\ 

ei L €l = e\l + z 1 fMU{e 1 ) dM {^) 



dN\ dN{e u M) 



deidt 



pz h (M,zi) 
X / dz" 

Jz^O 



dt 



dz" 



i,{z") . (14) 



Here z^,[M^Z\) is the maximum redshift a star may born, within comoving Hubble time 
1/H Q (1 + zi), to contribute to the local power depending on its mass. This can be calculated 
from equation ( flOl) by switching the sign of the second term inside the square brackets. 

Note that an integration of the luminosity density, the comoving quantity in equa- 
tion (Till) divided by (1 + zi) 4 and multiplying with it after integration, over cosmic time also 
leads to the comoving EBL energy density, an approach advocated by Salamon & Stecker 
(1998) and by Dwek et al. (1998). This amounts to adding up luminosity densities from 
all past epochs measured at a particular epoch to calculate the EBL. Both this approach 
and our formalism given by equations f|T2|) & ffT3|) produce the same EBL result as we have 
tested at different redshifts. 



3. Models of Star formation and dust attenuation 

The star formation history of the universe and the number of stars formed at a given 
epoch are not independent. A combination of both are needed to fit the SFR data (see, e.g., 
Madau et al. 1996; Cole et al. 2001; Hopkins & Beacom 2006) and the luminosity density 
data from local galaxy surveys (see, e.g., Baldry & Glazebrook 2003; Driver et al. 2008). 
The classic Salpeter (1955) IMF is still preferred by astronomers with some modifications. 
The modified "Salpeter A" model IMF is dN/dM ex M~ K with k = 1.5 below O.5M and 
2.35 above O.5M . Models by Scalo (1986; 1998) are in violation with a strong upper limit 
of k < 2.7 above 1M as found by Baldry & Glazebrook (2003). We also use the "Baldry- 
Glazebrook" IMF model with k = 1.5 below 0.5M© and k = 2.2 above O.5M . 

Cole et al. (2001) found a parametric form of the SFR given by ip{z) = h(a + bz)/[l + 
{z/c)% with parameters (a, b, c, d) = (0.0166, 0.1848, 1.9474, 2.6316). Hopkins & Beacom 



(2006) also used this parametric form, as well as their own piecewise fit in the form ~ 
10 a '(l + zf\ along with the "Salpeter A" and "Baldry-Glazebrook" IMF to fit SFR data. 
The best-fit parameters they found are given in their Tables 1 and 2. We define five models 
of different IMF and SFR combinations as: 

• Model A: Cole et al. (2001) SFR and Salpeter A IMF 

• Model B: Cole et al. (2001) SFR formula fitted by Hopkins & Beacom (2006) with 
Salpeter A IMF 

• Model C: Cole et al. (2001) SFR formula fitted by Hopkins & Beacom (2006) with 
Baldry-Glazebrook IMF 

• Model D: Hopkins & Beacom (2006) SFR with Salpeter A IMF 

• Model E: Hopkins & Beacom (2006) SFR with Baldry-Glazebrook IMF 

We plot these models in Fig. [31 Note that there is a significant amount of uncertainty 
among the models even at z ~ 0. Also, at the highest redshifts the SFR models may 
be underestimating the true rate (see, e.g., Faucher-Giguere et al. 2008). Most of the 
uncertainty comes from the dust in the host galaxies that absorbs stellar emission and reemits 
it into infrared wavebands, forming the low energy, < 0.1 eV, part of the EBL. A precise 
model of the fraction of starlight which directly escapes the host galaxy is still missing. Such 
a model would depend on the galaxy types and their orientations as well as on redshift. 

Bearing in mind the uncertainties in dust absorption models discussed above, we adopt 
a model recently developed by Driver et al. (2008), who have calculated the averaged photon 
escape fraction / esc (A) from observations of 10,000 nearby galaxies convolved with galactic 
dust models. We fit their results with four segments as 

{0.688 + 0.556 logA ; A < 0.165 
0.151 - 0.136 logA ; 0.165 < A < 0.22 
1.0 + 1.148 logA ; 0.22 < A < 0.422 ( ' 

0.728 + 0.422 logA ; A > 0.422 , 

with A in microns, and assume it to be universal or independent of redshift. The escape 
fraction / esc p» above ~ 10 eV as photons above this energy are absorbed by galactic 
gas. A feature of this model (Driver et al. 2008) is that the total amount of stellar energy 
absorbed by dust, in our local universe, is equivalent to the observed total luminosity density 
in infrared photons. The emission from luminous infrared galaxies, some of which may have 
an AGN core, is an added correction to the infrared EBL (see also Kneiske et al. 2004). 
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Results 



We use Models A - E for the SFR+IMF combinations and equations f fT2|) and (fT4j) . 
respectively, to calculate the EBL at z = and luminosity density at z — 0.1 in starlight 
component, and compare with data. The results, with numerical integrations carried out by 
the multidimensional adaptive Monte Carlo code Vegas (Lepage 1978; 1980), are plotted in 
Figs. H] and [5] for the SPL (left panels) and BPL (right panels) models of the L — M relation. 

To plot the luminosity density data points, in Fig. HI we took the AB magnitude data 
rriAB = j + 2.5 log h with j = —2.5 log(L e /W Hz~ 1 Mpc -3 ) + 34.1 and converted them to 
luminosity density as eL t = /i(c/A)10 13 - 64_mAS / 2 - 5 W Mpc~ 3 . The factor 34.1 comes from 
2.5 log(4.345 x 10 13 ), where 4.345 x 10 13 W Hz" 1 Mpc" 3 is the minimum of the absolute AB 
magnitude scale. Note that all data points are not measured at the same redshift or corrected 
to z = as they come from different surveys: The Sloan Digital Sky Survey (SDSS), the Two 
Micron All Sky Survey (2MASS), the Two-Degree Field Galaxy Redshift Survey (2dFGRS) 
and the Galaxy Evolution Explorer (GALEX). The data points by Budavari et al. (2005) 
are calculated for redshift 0.07 - 0.13 and the mean redshift is 0.1 for the data points by 
Blanton et al. (2003). The data point by Norberg et al. (2002) and those by Cole et al. 
(2001) are corrected to z = 0. The redshift for Kochanek et al. (2001) data point is 0.03. 
Note that we plotted our luminosity density models A-E for redshift 0.1 to be consistent 
with most data points. 

The three EBL data points in Fig. [5], in the UV band, are calculated by Bernstein, 
Freedman & Madore (2002) using measurements from the Hubble Space Telescope (HST). 
However, revised (Mattila 2003; Bernstein, Freedman & Madore 2005) and re-revised (Bern- 
stein 2007) estimates have put substantial uncertainty on these data points. The rest of 
the data points come from analyzing measurements by the Diffuse Infrared Background 
Experiment (DIRBE) onboard the Cosmic Background Explorer (COBE) satellite by differ- 
ent authors (Dwek & Arendt 1998a; Gorjian, Wright & Chary 2000; Wright & Reese 2000; 
Cambresy et al. 2001 and Levenson, Wright & Johnson 2007). This is the reason for multiple 
data points at the same energies. Note that the ~ 1 eV (1.25 /im) data point by Levenson, 
Wright & Johnson (2007) is really a la limit. We plot only those upper limits (Dwek & 
Arendt 1998a; Hauser et al. 1998) and lower limits (Madau & Pozzetti 2000; Fazio et al. 
2004) which are directly derived from the DIRBE and HST data. 

Given a combination of mixed data points, Models A, B and D fare well reproducing 
the luminosity density data in Fig. H] with the SPL model of L — M relation (left panel) given 
by equation (jl]). These models miss the two data points at the smallest and largest (within 
2a) wavelengths, however. In case of the BPL model of L — M relation (right panel) given 
by equation §5$), Models A, B and D agree with two small wavelength data points better 
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than the SPL model, but miss most data points at longer wavelengths. Models C and E are 
systematically lower than all data points for both the SPL and BPL cases. The integrated 
total energy output by stars, using Model B (SPL), is 1.2 x 10 35 W Mpc -3 before dust 
absorption and 0.7 x 10 35 W Mpc -3 after dust absorption. As for comparison, Driver et al. 
(2008) calculated these values as 1.6 x 10 35 W Mpc -3 and 0.9 x 10 35 W Mpc -3 , respectively. 

Comparing with EBL data in Fig. [5j we find again that Models A, B and D with the 
SPL L — M relation (left panel) represent the data points better than all other models in 
both SPL and BPL cases. Models A, B and D in the BPL case are consitently lower, below 
2 eV, than the SPL case, a reflection of the trend observed in the luminosity density plots 
as well. Models C and E in both the SPL and BPL cases (right panel) are below the HST 
lower limits (Madau & Pozzetti 2000). None of our models are able to reproduce EBL data 
below 1 eV hinting that an extra component is required for modeling infrared data. 

Models A, B and D involve Salpeter A IMF which predicts more stars forming below 
~ 2M Q than the Baldry & Glazebrook (2003) IMF used in Models C and E. This may 
be one reason why these models underproduce the local luminsoity density and EBL, since 
stars with < 2M dominantly produce < 2 eV photons. The other affecting factors are the 
differences in SFR models (see Fig. [3]) and mass-to-light ratio. With fi > 1 in equation (jSJ), 
for the Baldry-Glazebrook IMF, it may be possible to close gaps between Models C and E 
with other models with BPL L — M relation. However, we do not explore that possibility in 
the present paper. For the purposes of 7-ray astronomy, any of the Models A, B or D can 
reproduce local UV/optical EBL data. Next we compare our EBL models to the models by 
other authors. 



4.1. Comparisons with other authors 

Among several existing models of the EBL, the one by Primack, Bullock & Somerville 
(2005), plotted in Fig. [6] is used to be considered as the "Low" EBL model. It is consistent 
with the lower limits from the galaxy counts above ~0.5 eV. An updated version of the 
fast evolution model by Stecker, Malkan & Scully (2006) is also plotted in Fig. [6j This 
is considered to be the "High" EBL model. The best-fit model by Kneiske et al. (2004), 
plotted here, and the phenomenological fits by Dermer (2007), not plotted here, are generally 
in between these two models. 

Primack, Bullock & Somerville (2005) used Monte Carlo simulations of galaxy evolu- 
tion and emission by the evolving galaxy population to calculate their EBL model. Kneiske, 
Mannheim & Hartmann (2002); and Kneiske et al. (2004) used the outcome of the simu- 
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lations of stellar population luminosity at different redshift, and integrated over redshift to 
calculate the EBL. Stecker, Malkan & Scully (2006) assumed that the observed luminosity of 
a galaxy at 60 fim can be used to calculate its luminosity at all wavelengths, and used galaxy 
luminosity functions to calculate the EBL. In our models, we used blackbody emission from 
stars, based on the solar models, and convolve with the star formation and initial mass func- 
tion models to calculate the EBL. These methods, therefore, are not directly comparable to 
each other. However, assumptions about an SFR or/and an IMF model(s) should affect all 
methods as also shown by Kneiske, Mannheim & Hartmann (2002); Kneiske et al. (2004). 
Indeed the variations between our EBL models A-E results from the differences in SFR and 
IMF combinations. 

Our models (A, SPL), (B, SPL) and (D, SPL) are consistent with the HST lower limits 
above ~1 eV and agree with the EBL model by Primack, Bullock & Somerville (2005) in the 
~l-3 eV range. At energies lower than ~1 eV our models are lower than all other models. 
This is probably because the infrared EBL component, from dust radiation; luminous in- 
frared galaxies and post main-sequence stars, becomes important at these lower energies. At 
energies above ~3 eV all other EBL models are higher than our models. There may be ad- 
ditional contributions from AGN and white dwarfs at these higher energies, however due to 
a lack of data points comparisons between the models above ~3 eV become less meaningful. 

Both the models (B, SPL) and (D, SPL) better represent the luminosity density data in 
Fig. IHthan the model (A, SPL) and can be used as our best-fit EBL models. For illustration 
and further calculations, we have chosen to use Model B, SPL plotted in Fig. We have 
also plotted our Model (B, BPL) for comparisions. 



4.2. Evolution with redshift 

Understanding the evolution of the background light with redshift is a key to 7-ray 
astronomy, which we discuss in the next section. The mean-free-path for 77 absorption with 
the EBL spans astronomical distances, over which the EBL itself changes noticeably. We 
have plotted the comoving EBL energy density at different redshifts 0-5 in Fig. [7] using equa- 
tion (Tl3|) for the Model B, SPL. Initially the EBL density increases with redshift because of 
the sharp rise in star formation (see Fig. [3]) below z ~ 2 and a decreasing volume. At redshift 
> 2, the EBL density decreases because the total number of stars formed up to that redshift 
from z = 6 decreases. The contribution to the EBL above ~2-3 eV dominantly comes from 
high mass stars. Since the livetimes of high-mass stars are shorter, their contribution to 
the EBL at a particular z is mostly determined by how many of them are formed at ~ z 
(see Fig. [T]). On the other hand the overall population of low mass stars, which dominantly 
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produce lower energy photons, increases with decreasing redshift. As a result, the ratio of 
high-energy photon density increases with redshift (as evidenced by the hump at ~7 eV) 
compared to the low energy photon density reflecting a decreasing overall population of low 
mass stars for z — » 6. 



5. Implications for 7-ray Astronomy 

High-energy photons in the ~10-300 GeV energy range from sources such as GRB at 
* > 0.5 are subject to 77 — > e + e~ absorption by the EBL starlight photons as we modeled 
here. To calculate the 77 absorption opacity one needs to take into account an evolving 
EBL with redshift as plotted in Fig. [7J We calculate this opacity both numerically, using 
interpolation to the exact results from Model B, SPL calculation, and analytically, by fitting 
the Model B, SPL results. 

We provide such a fit below and plot it in Fig. [S] at different redshift overlaid with 
numerical calculation. The polynomial fit parameters are 

log(eiit £1 /ergs cm" 3 ) = A + A x x + A 2 x 2 + A 3 x 3 + A 4 x 4 

x = log( ei (eV)/(l + *i)) 
A Q = -14.4829 + 0.8275zi - 0.2451*? + 0.0046* 3 + 0.00434 - 0.0004?? 
A 1 = 0.3157 - 1.105,2! + 1.1026*? - 0.4764* 3 + 0.092? - 0.0062*? 
f -1.9888 + 1.6527*1 + 1.0294*? ; *i < 0.8 
2 ~ \ -0.5549 - 0.0295*i - 0.1133*? + 0.0079z 3 ; z x > 0.8 
A 3 = -0.1507 + 0.9114*! - 1.8907*? + 0.8816* 3 - 0.1837*? + 0.0141*? 
A 4 = 0.3014 - 3.5371zi + 2.4574z? - 0.8474* 3 + 0.1343*? - 0.0079*? , (16) 

where e\ = e(l + *i). The fit is generally good at the tens of percentage level of the numerical 
results for e\ > 0.5 eV for all redshift z\ < 5. 

With a fit to the diffuse background photons, in equation f[T6|) . we can easily calculate 
the optical depth of 77 absorption for an energetic photon originating at redshift * with 
observed energy E as (Gould & Schreder 1967; Brown, Mikaelian & Gould 1973) 

dt 



r 77 (E,*) = c / dz\ 
Jo 



dz\ 

x / dt\ I dcosO — — (1 — cos6')cr 77 (s) 
'0 J-i 2 €i 



cirr„ 



m^c 8 f z dz\ 



E 2 Jo (1 + ziY 



dt 



dz\ 



13 



de^Mei)]- (17) 

I mlc A / E(l+ Zl ) e l 

for an isotropic background photon field. Here cr 77 (s) is the total 77 — > e + e~ cross-section 
(see, e.g., Jauch & Rohrlich 1955) and s = E{1 + Z\)€\{1 — cos6')/2mgC 4 is the center-of- 
mass energy squared. The function <p[so], with so = E(l + z^ei/m^c 4 , is given in Gould & 
Schreder (1967). The threshold energy for e + e~ pair production from the condition Sq = 1 
is ei jt h ~ m 2 e c A / E(l + z\). The lower limit of the energy integration ei^h ~ 1 eV and 
~ 250/(1 + z\) GeV 7-rays interact dominantly with these photons. Higher energy 7-rays 
mainly interact with softer photons, at an energy range where the re-processed dust emission 
dominates the EBL and our models do not fit the data. With an upper limit ei jmax ~ 10 eV, 
it is safe to use equation ( TT7I) in the ~10-300 GeV 7-ray energy range for redshift > 0.5. Note 
that this the relevant energy range for the Large Area Telescope on board the Fermi Gamma 
Ray Space Telescope. Air Cherenkov Telescopes also become sensitive to 7-rays >50 GeV. 
The survival probability for a 7-ray created at redshift z to reach Earth is exp[— r J7 (E, z)]. 
The observed flux of 7-rays is attenuated by this factor from the source flux. 

We have plotted the 77 opacity (r 77 ) in Fig. [9] for redshift 0.5-5.0 and observed 7-ray 
energy 10-300 GeV using equation ( fT7j) . The solid and dashed lines correspond to the exact 
calculation and calculation using the fit in equation ([TBI) respectively. As can be seen, the 
two methods of calculation give results within tens of percent. For z > 3 (right panel) 
we have divided r 77 with constant factors given on the plot to avoid cluttering of different 
curves. We also provide our exact calculation of r 77 in Table [JJ for different redshift and 
7-ray energies. The relation r 77 = 1 (Fig. [TO]) , sometimes referred to as the Fazio-Stecker 
relation (Fazio & Stecker 1970), corresponds to a 7-ray horizon of the universe. 

With the calculated opacities, we can calculate the true 7-ray flux of a distant source 
before absorption in the EBL by multiplying the observed flux /# with the factor exp[r 77 (£ 1 )]. 
In case of the furthest known blazar 3C279 at z = 0.536, the flux data points measured by 
MAGIC (Albert et al. 2008) are plotted in Fig. [UJ We have also plotted the deabsorbed 
data points (open circles) using our EBL Model B (SPL) and a power-law fit dN/dE = 
N (E/200 GeV)- r with N = 1.5 x 10~ 9 TeV" 1 cm~ 2 s" 1 and T = 2.8. The deabsorbed 
spectrum of 3C 279 thus becomes harder than the observed value of T = 4.11 (Albert et al. 
2008). Note, however, that our 77 opacity calculation is less reliable for the highest energy 
data point at 474 GeV. The T for the deabsorbed spectrum, thus, can be < 2.78 if r 77 is 
much larger than our estimate. 

Simple Fermi mechanisms of particle acceleration in both relativitic and non-relativistic 

shocks, result in particle spectra dN/dE oc E~ p with p w 2. In the case of a one-zone 
synchro-Compton mechanism to produce TeV 7-rays at the shocks, the source spectrum of 
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7-rays typically should be larger than T = (p + l)/2 ~ 1.5. Spectra harder than the simple- 
minded limiting value may arise in a number of scenarios. For a special environment at or 
surrounding the shocks, effects such as a cutoff of shock-accelerated elecron spectrum below 
a very high energy ~GeV (Katarzyhski et al. 2006; Stern & Poutanen 2008) or internal 77 
absorption in a dense region of quasi mono-energetic soft photons (Aharonian, Khangulyan 
& Costamante 2008) may result in a very hard spectrum T < 1.5 over a short energy range. 
Acceleration at shear flows (Stawarz & Ostrowski 2002; Rieger & Duffy 2006) or Monte Carlo 
shock-acceleration models may also result in T < 1.5 (Virtanen & Vainio 2005; Stecker et al. 
2007). 

Models with multiple spectral components, arising from multiple zones, are attractive 
to explain observed hard blazar spectra. Recently Bottcher, Dermer & Finke (2008) have 
invoked such a model of Compton upscattering of CMB photons in large scale jet to explain 
hard spectrum of 1ES 1101-232. Also, cosmic-ray acceleration and interactions at the shocks 
and subsequent cascade in the EBL-CMB may also produce very hard 7-ray spectra (Coppi 
& Aharonian 1997). However, such models cannot explain highly variable TeV spectra, e.g., 
from 3C 279 (Albert et al. 2008). 



6. Conclusions 

We have derived a class of well-defined models for the spectral energy density of the 
EBL. The models with modified Salpeter A initial mass function, a single power-law mass- 
luminosity relation and Cole et al. (2001) or Hopkins & Beacom (2006) star formation history 
reasonably fit the EBL UV-optical data and the luminosity density in our local universe. 

Our models are based on the underlying assumption that the bulk of the EBL radiation 
between ~ 0.1 - 10 eV is due to stellar radiations absorbed by dust, which can be determined 
from recent analyses of galaxies by Driver et al. (2008). This approach differs from models 
by Stecker, Malkan & Scully (2006) based on luminosity evolution of galactic spectral energy 
distributions, which is limited to the accuracy of the available observational data used in the 
survey. Without need of a population synthesis code (e.g., models by Primack, Bullock & 
Somerville 2005) or fits to the results of such a code (e.g., models by Kneiske, Mannheim & 
Hartmann 2002), our model is based on well-studied results from stellar astronomy. 

The sources of uncertainties in our models are the (i) main-sequence age of the star 
and stellar luminosity, discussed earlier, (ii) star formation history and initial mass function, 
and (iii) dust absorption. We have already discussed point (ii) in some details using five 
SFR+IMF models. Note that in all those models, the IMF was assumed to be independent 
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of redshift. In principle the normalization of the IMF or even its shape may depend on 
z. Nevertheless, a universal IMF fits SFR data reasonably well. The evolution of the dust 
absorption model with redshift, which we have not taken into account, is of potentially 
greater concern. At high redshift it is more reasonable to assume that the dust absorption 
(escape) fraction would be higher (lower), so that the stellar contribution to the EBL at high 
z would be less than if using a constant absorption fraction. A more detailed examination 
of these issues are under further study (Finke et al., in perperation). 

We have provided an analytic fit to our best-fit EBL model and its evolution with 
redshift. This result can be used to calculate, as we have done in this work, opacity of the 
universe to ~10-300 GeV 7-rays relevant for the high energy data from the Fermi Gamma 
Ray Space Telescope and Air Cherenkov Telescopes, and estimating unknown quantities such 
as the spectrum and energy at production of the distant GRBs and blazars such as 3C 279. 

We thank Eli Dwek, Claude-Andre Faucher-Giguere, Dieter Hartmann, Tanja Kneiske, 
Kalevi Mattila and Floyd Stecker for helpful comments. The work of S.R. and J.D.F. was 
supported by the National Research Council Research associateship program at the Naval 
Research Laboratory. The work of C.D.D. was supported by the Office of Naval Research. 
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Table 1. 7 ray absorption opacity (r 77 ) with EBL Model B, SPL at different redshift 
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Fig. 1. — Plot of the redshift (max{0, z<i(M,z)}) at which a star of mass M which was 
born at redshift z =1 - 6 evolves off the main sequence using equation ffTUl) and assuming 
the standard (0.7, 0.3, 0.7) ACDM cosmology. The two sets of curves correspond to SPL 
(dashed) and BPL (dotted) L — M relation, respectively, in equations (j4j) and (JSJ). 
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Fig. 2. — Initial mass function (IMF) models, assumed to be universal. The integral of 
dN/dhiM is set to unity over M=(O.1-12O)M range. 
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Fig. 3. — Models of SFR combined with specific IMF. Model A correspond to Cole et al. 
(2001) SFR with Salpeter A IMF, Model B correspond to Cole et al. (2001) SFR formula 
fitted by Hopkins Sz Beacom (2006) with Salpeter A IMF, Model C correspond to Cole et 
al. (2001) SFR formula fitted by Hopkins & Beacom (2006) with Baldry-Glazebrook IMF, 
Model D correspond to Hopkins & Beacom (2006) SFR with Salpeter A IMF and Model E 
correspond to Hopkins & Beacom (2006) SFR with Salpeter A IMF. 
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Fig. 4. — Total energy output or luminosity density, e\L ei from equation (fi~4|) . in direct 
starlight by our local (zi = 0.1) universe. The data points are from SDSS, 2MASS, 2dFGRS 
and GALEX surveys of local galaxies. The smooth lines are our calculations using equa- 
tion (fi~4|) and based on SFR+IMF models A-E described in Section |3j The left and right 
panels correspond to the SPL and BPL models of L — M relation, respectively, given by 
equations (T4J) and ([5]). 



Wavelength (u.m) 
1 




Wavelength (flm) 
1 




Energy (eV) 



Energy (eV) 



Fig. 5. — EBL energy density, eu e = e 2 dN(e,z = 0)/dedV from equation ( TT2|) . from direct 
starlight in our local universe. The smooth lines are calculated using equation (fT2|) and 
correspond to SFR+IMF models A-E as defined in Section[3l The data points with errorbars, 
and lower (triangles) and upper (inverted triangles) limits are described in Sec. HI Left 
panel — models with L — M relation (SPL) in equation (jlj). Right panel — models with 
L — M relation (BPL) in equation ([5]). 
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Fig. 6. — A comparison of the EBL models (B, SPL) and (B, BPL) derived in this paper with 
models by other authors. The updated fast evolution model by Stecker, Malkan & Scully 
(2006) is based on backward-evolution models of local galaxies. The models by Primack, 
Bullock & Somerville (2005) is based on Monte Carlo simulatons of galaxy evolution with 
initial conditions. The model by Kneiske et al. (2004) is based on results of population 
synthesis models. Our models are based on initial mass function and star formation models 
convolved with stellar properties. 
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Fig. 7. — Evolution of the comoving EBL energy density e\U ei as a function of the comoving 
photon energy t\ with redshift for our best-fit EBL model (Model B with SPL L — M 
relation). Model (D, SPL) gives simililar results. The curves are plotted for z = 0-5 with 
0.2 interval. The density increases first and then decreases as indicated by the solid and 
dashed arrows. 



-27- 




Energy (eV) 



Fig. 8. — Fits (smooth lines) to an evolving EBL energy density e\U ei (points) as a function 
of the comoving frame photon energy ei at different redshifts for Model (B, SPL). The fit 
function is given in equation ( fT6l) . 
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Fig. 9. — Gamma-ray opacity of the universe at different redshifts (Left panel: 0.5-2.5, Right 
panel: 3.0-5.0) and energies as plotted here. The solid and dashed lines correspond to the 
exact calculation and calculation using the fit to the EBL in equation [16], respectively. We 
rescaled r 77 in the Right panel to separate diferent curves. The numerical values for the 
exact calculation are also listed in Table [IJ 
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Fig. 10. — Gamma ray opacity contrours in the E — z plane using the EBL fit in equa- 
tion (TT61) . The r 77 = 1 contour plotted here is known as the Fazio-Stecker relation (Fazio &: 
Stecker 1970) and represents a 7-ray horizon of the universe. 
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Fig. 11.— High-energy 7-ray flux from the furthest known (z = 0.536) blazar 3C 279 
as measured by MAGIC (filled circles with errorbars. We use our EBL Model B, SPL to 
estimate (deabsorb) the source flux (open circles with errorbars) by multiplying the observed 
data points with exp(r 77 ). The dashed line is a power- law fit dN/dE ex E~ v to the deabsrbed 
data points resulting in a spectral index V = 2.8. 



